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Abstract 

In this paper, the scattering/transmission inside a step-modulated subwavelength metal slit is 
investigated in detail. We firstly investigate the scattering in a junction structure by two types of 
structural changes. The variation of transmission and reflection coefficients depending on structural 
parameters are analyzed. Then a multi-mode multi-reflection model based on ray theory is proposed 
to illustrate the transmission in the step-modulated slit explicitly. The key parts of this model 
are the multi-mode excitation and the superposition procedure of the scatterings from all possible 
modes, which represent the interference and energy transfer happened at interfaces. The method we 
use is an improved modal expansion method (MEM) , which is a more practical and efficient version 
compared with the previous one [Opt. Express 19, 10073 (2011)]. In addition, some commonly 
used methods, FDTD, scattering matrix method, and improved characteristic impedance method, 
are compared with MEM to highlight the preciseness of these methods. 

PACS numbers: 78.68. +m, 78.20.-e, 42.79.Gn, 73.20.Mf 
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I. INTRODUCTION 



Subwavelength metal slits, as a kind of metal/insulator/metal waveguides, have attracted 
much attention in recent years not only because of their ability to guide light beyond the 
diffraction limit, but also because of several remarkable advantages, such as strong field 
localization, simplicity, and convenience for fabrication and integration into optical circuits 
[1-14]. When light (infrared, visible spectrum) propagates along a metal/air interface, it 
will excite a collective oscillation of free electrons at the surface of the metal, causing a 
field exponentially decaying away from the interface. This mode is called as surface plasmon 
polariton (SPP) [2-6]. In a subwavelength metal/air/metal slit, the case is somehow different, 
since the in-slit SPP [5] wave decays exponentially in the metals and is fiat in the air, which 
is the lowest eigenmode in the sht structure and the core part of the subwavelength metafile 
optics. 

Step modulation is one of the key elements in photonic engineering that are employed in 
subwavelength metal structures to design and fabricate functional plasmonic devices, such 
as filters [7-10], refiectors [11], and photonic bandgap structures [11,12]. Besides, the step 
modulation is of important theoretical significance since they are helpful for investigating 
SPP scattering. Had the knowledge and combined with the staircase approximation and 
transfer matrix technique, numerical results of more complicated structures can be obtained. 

Up to now, a number of methods have been used to calculate the SPP scatter- 
ing/transmission inside a step-modulated slit. The finite-difference time-domain method 
(FDTD) is a well developed simulation method that provides relatively accurate results and 
has been considered as a standard for testing other theoretical methods [5-14]. The effective 
index method, on the other hand, is a simplified and direct theoretical method where only 
the SPP modes are involved in the calculation, a method of one mode approximation; how- 
ever, this simplification causes the loss of the scattering information considerably and the 
numerical imprecision turned out to be considerable under some conditions [6]. Matsuzaki 
et al. presented a transmission model and gave a better description of the SPP scattering 
using the characteristic impedance method [13]. Pannipitiya et al. [14] suggested an im- 
proved version of this method, which will be called as improved characteristic impedance 
method (ICIM) in this paper. Lin et al. presented a similar transmission model and used 
the scattering matrix method (SMM) to calculate the transmission [7-9]. Although the cal- 
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culated results from these two methods fit the FDTD results, two approximations, one mode 
approximation and quasi-statistic approximation, are used in the calculation, which limits 
their application scope and numerical precision, as will be discussed in Sec. 3 below. Re- 
cently, we successfully applied the modal expansion method (MEM) in discussing the wave 
behavior inside a symmetric step-modulated slit [6]. This method did not involve the two 
mentioned approximations and provided more accurate results. 

In this paper, the MEM is further improved so as to apply to the asymmetric modulated 
case in investigating the scattering/transmission mechanisms inside a step-modulated slit. 
A multi-mode multi-reflection model is proposed to explain the transmission process. A 
remarkable advantage of MEM is that its precision is controllable. This enables us to discuss 
the preciseness of FDTD and ICIM by comparing the results from these methods and MEM. 

The paper is arranged as follows. Section 2 sets our model of a step-modulated metal 
slit and presents the improved MEM formulas. In Sec. 3, the scattering in a junction 
structure is studied firstly as a prerequisite for later discussion, and then a multi-mode 
multi-refiection model is proposed to reveal the transmission mechanism in a step-modulated 
slit. Comparisons between MEM, FDTD, SMM, and ICIM are also given in this section to 
highlight the restriction of the one mode approximation and quasi-statistic approximation. 
Finally, conclusions are presented in Sec. 4. 

II. MODEL AND THE IMPROVED MEM 

In this section, we set the model of single-sht structure and present the formulas of the 
improved MEM which is more practical and efficient compared with that in our previous 
work [6]. 

The model structure is shown in Fig. 1. It is infinitely large in yz plane but confined 
in X direction by perfectly conducting walls at x = and x — L. The structure is divided 
into three regions along x direction, composed of silver /air /silver, and three layers along 
y direction. The lower boundary of the Zth layer is labeled as Q^'"^^ and the interfaces 
between regions are by xf^ and x'2 ■ The denotations g^'^ and w^^'> label layer height and 
slit width, respectively. A transverse magnetic (TM) wave, with magnetic field H being in 
z direction, is normally launched from y = Q^'^^ in Layer 1 and propagates upward. The 
structural parameters are given in the caption of Fig. 1. 
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FIG. 1: Sketch of a step-modulated metal slit structure confined in x direction with perfectly 
conducting walls at and L = 2 /xm. A TM wave with wavelength Aq is normally launched at 
y = Q(0). Q(l)=0. 



The dielectric constant of silver as a function of the wavelength of the incident wave Aq 
is evaluated as e^g = (3.57 — 54.33Ao) + 0.083Ao + 0.921Aq) by fitting the experimental 
data [15], which is valid for 0.6 < Aq < 1.6 /xm. In this paper, the wavelength is mainly set 
as Aq = 1 /jLia; thus, EAg — —50.76 + 0.083i. 

In the remaining part of this section, we suggest an improved version of the MEM, which 
has the same output as the previous one [6] but is easier and faster. 

The substance of MEM is to expand the unknown functions (electromagnetic field distri- 
bution in present case) by a complete set of orthogonal functions. This makes the MEM has 
two folds: one is the eigenvalue problem of the system; the other is to establish and solve the 
coupled equations subject to boundary conditions. However, the choice of the complete set 
in modal expansion is not unique, but depends on the configuration of the system. It can 
be the eigenfunctions of a specific structure or other functions such as sine or exponential 
functions. 

In Ref. [6], the fields in the given structure were handled by separation of variables. 
The factors containing x variable were expanded by the eigenmodes {'4'n\^)} between two 
perfectly conducting walls. The magnetic fields were expressed as 



g(°) <y< g(^) 

g(i) <y< Q(2) 
g(2) < y < OO, 



(1) 



where /„, Rn, En, Fn, and Tn were expansion coefficients which involved the scatter- 
ing/transmission information of every eigenmode. It was inevitable to solve a transcendental 
equation in order to achieve the eigenvalues and eigenf unctions. Even with assistance of a 
powerful root-seeking method [16], this procedure was still time-consuming. Moreover, each 
layer had its own eigenfunctions. At an interface, it was required by the boundary conditions 
to calculate the overlap between the eigenfunctions at the two sides of the interface, called 
as coupling integrals. Such integrals brought complexity to the program. 

To avoid these difficulties, in this paper the factors containing x variable were expanded 
by a sine basis subject to the perfectly conducting boundary condition. That is to say, the 
complete set {ipn{x)} is chosen as 



(2) 



The eigenvalues kxn are solely determined by the distance between the two perfectly con- 
ducting walls, independent of the positions xf^ and X2\ so that is valid for all the three 
layers. 

Correspondingly, the magnetic fields and its derivative in the three layers can be expressed 
as [17] 



Hz{x,y) = < 



(3), 



<y< g(i) 
g(i) <y< g(2) 

g(2) < y < oo. 



(3) 

where im, r^, e^, fm, and tm are the expansion coefficients. The insertions of Eq. (3) 
into Helmholtz equation yields an eigenvalue problem in each layer expressed by ^(')iy(') ~ 
{ik^y ^y^W^''^ with ik^^^ and W^^^ being the eigenvalues and eigenfunctions, and the operator 
being [18] 

A« = - |a;2[/] + [K] [i^]| , (4) 

where 



E^''^] =/o^v'„(x)^„(a;)£«(a;)dx 

J mn 

[K]mn = lo ^mix)-^(Pn{x)dx 
[-^Imn ~ ^mn- 



(5) 
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In these equations , [•] denotes a N x N matrix where N is the truncation number, and 
ko — 2t:/\q is the wave vector in vacuum. 

Here we mention the two advantages of the sine expansion in x direction. One is that the 
eigenvalue problem, Eq. (4), is very easy for computer implementation, which avoids the 
cumbersome solution-seeking procedure necessary in the eigenmode expansion [6]. The other 
is that the complete set of sine functions are the same for all layers, so that the coupling 
integrals at the interfaces become quite simple. These two advantages make the calculation 
program greatly simplified. 

Corresponding to Eq. (3), the derivative of the magnetic field is expressed as 



1 d 



e dy 



Hz{x, y) 



M^y- Ti^(2) •7,(2) 
Z-^n £(2) Z^m " n,m''^ym 



Cm ,6 



(3) , 



Q(o) <y< g(i) 

g(^) <y< Q(^) 

< y < oo. 



(6) 



Applying the layer boundary conditions, we obtain the coupled equations as follows: 



i—im ' ' pni 



pn 



Z^m ' ' pm 



pn nm ym 



m ' pm 

y E^^^y py(2)^(2) 

pn ' ' nm ym 



+ fr, 



J m^ 

^y ^(2)^^(2) 



i-^m ' ' pm 



(7) 



,1.(2) „(2) „ 



V M/(3)A;(3)t 



nm ym 



The incident coefficients i^n are determined by the incident wave. In this paper, the 
incident wave is always a SPP wave launched in Layer 1, namely, '^^^{x). Therefore, one 
naturally has 

[%p{x)4'\x)dx = E^^m^^m, (8) 
■^0 m 

which determines the coefficients im- After setting the incident coefficients i^, the four 
groups of coefficients, r^, e^, fm, and t^, can be obtain from Eq. (7). Thus all the field 
quantities are obtained. 

In this paper, we will focus on discussing the reflection/transmission mechanisms, which 
arc mainly presented by the reflection coefficients and transmission coefficients T„ in 
Eq. (1). For example, the amplitudes of the SPP modes in Layer 1 and 3, and |T„|, 
are the reflection and transmission efficiencies of the system, and their arguments, arg{Rn) 
and arg{Tn), are the corresponding phase shifts. Therefore, a projection between the fields 
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calculated by Eq. (3) and the eigenmodes {i/^^^x)} is implemented for obtaining the i?„ 
and T„. In the following, the absolute values of these coefficients may generally be named 
as excitation efficiency. 

Thus we accomplish our formulation presentation. This method is briefly outlined as 
follows: the flcld is expanded by sine functions which are complete and uniform in all layers. 
The corresponding eigenvalue problem becomes a matrix form as shown in Eq. (4), which 
makes the calculation quite easy. Accordingly, the procedure here is much more practical 
and efficient compared to the previous one [6]. At last, the overlaps between the calculated 
field and the eigenmodes in Eq. (1) give the required refiection and transmission coefficients 
necessary for physical analysis. 

Although we merely study the three-layer structure, our procedure developed here is 
easily applied to more complicated structures by implanting the S matrix algorithm [19] or 
the enhanced transmittance matrix approach [20]. 

III. NUMERICAL RESULTS AND ANALYSIS 

In this section, we investigate the scattering/transmission mechanisms inside a step- 
modulated subwavelength metal slit. To do so, the scattering in a junction structure is 
first discussed in detail, since the slit comprises more than one junction structure. Then we 
disclose the multi-mode multi- refiection model in the transmission process in the slit. By 
the way, the numerical precision of FDTD and ICIM is discussed by comparing results of 
these two methods and MEM. 

In calculation, the confinement is set as L = 2 /xm. We have tested that 800 modes, 
= 800 in Eq. (4), are enough to give results with precision up to four significant digits. 
The convergence test for truncation number A^ and the preciseness test for confined width 
L will be carried out later in Fig. 7. 

A. Junction structures 

A junction structure is the connection of two half-infinitely long slits with widths being 
denoted by w^^^ and w^^\ respectively, which can be easily realized in Fig. 1 by setting the 
height of Layer 2 to 0. In Ref. [6], some scattering properties of symmetric structures have 
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been revealed. For example, the main ingredient of the fields inside the slits were guided 
modes which played a very important role in scattering/ transmission, and the unimportant 
components were the radiation modes excited which were necessary to fulfill boundary con- 
dition, but had little contribution to the transmission. So the following discussion will focus 
on the guided modes. However, the discussion in Rcf. [6] was limited to the symmetric case. 
Here we present a detailed investigation on how the scattering is affected by asymmetry. 

Two types of structural changes are considered. In Type I, the widths of the two slits are 
fixed and the position of the narrower one can be anywhere between left to right, as shown 
in the inner panel of Figs. 2(a) and (c). In Type II, the left walls of the slits are aligned and 
the width of the narrower one is fixed, but that of the wider one can vary, as schematically 
shown in the inner panel of Figs. 3(a) and (c). 

The results of Type I structure are plotted in Fig. 2. The left wall of the wider slit is at 
Xi — 0.6 //m. For w^^^ — 0.1 and w^^^ — 0.8 //m, when the position of narrower slit is moved 
from the left to right, the excitation efficiencies and their phase shifts are plotted in Figs. 
2(a) and (b), respectively. In Figs. 2(c) and (d) are the excitation efficiencies and their 
phase shifts of structure w^^-' = 0.8 and w^^^ = 0.1 fxm. In Figs. 2(e) and (f) are the absolute 
values of the first three cigcnf unctions of the narrower and wider slits, respectively. In Figs. 
2(a), (c), (e) and (f) the absolute values are plotted because these quantities are complex. 
The eigenmodes in the narrower and wider slits are denoted by ■0^""^ ^-nd respectively. 
In Figs. 2(e) and (f), the lowest modes ijji"'"'^ and plotted by the black curves, are just 

SPP modes, and the second modes ■02^"^ and ■02^"*^ plotted by the red curves, are of actually 
antisymmetric wave functions within the shts. Note that the curve of IV'a'"*'*! is divided into 
three parts by two zeros. The sign of the central part of is contrary to the other parts. 

We notice that under our present parameters, only the first two eigenmodes of the nar- 
rower slit -^i""-* and '^2""'* within the slit (guided modes), see the black and red lines 
in Fig. 2(e). The wave functions of the higher modes mainly distribute within the metal 
(radiation modes), see, as an example, the third mode ■03^"^ iii Fig- 2(e). The behavior of 
I '03""^ I with the position of the narrower slit has to be explicitly given as following. When 
the slit moves rightwards, the width of metal at the right side of the slit becomes thinner, 
so that the hill is compressed. If the narrower slit is on the right side of the center position 
X — 1 /im, the hill will appear at the left side of the slit. While if the slit is just at or very 
near the center x — 1 /im, there will be two hills at the two sides of the slit, respectively. 
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FIG. 2: (color online). Scattering in Type I junction structures. The left wall of the wider slit is 
fixed at xi = 0.6 fim. 

g(0) ^ Q{i) ^ g{2) ^ 0. In ah the fi eures, the results of the wider slit are 
plotted by solid lines and those of the narrower slit by dash-dotted lines. The narrower slit moves 
from the left to right. In Figs, (a) to (d), the x-axes are its central position. For the structure 
with w^^^ = 0.1 and w^^'^ = 0.8 fim, (a) excitation efficiency; (b) phase shift. For the structure with 
w (1) = 0.8 and w^'"^^ = 0.1 fim, (c) excitation efficiency; (d) phase shift. In (e) and (f) plotted are 
the absolute value of the eigen functions in the narrower and wider slits, respectively, when their 
left walls are aligned. 



since the structure in this case is symmetric [6]. 

Three obvious features of excitation efficiencies can be seen in Figs. 2(a) and (c). The 
first is that all the curves there exhibit a central symmetry, because all the configurations 
are symmetric with respect to the central line at x = 1 fim. The second is, from comparison 
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of the black and red solid lines in Figs. 2(a) and (c), that the excitation efficiencies of the 
SPP modes in narrow shts are much larger than those of the second modes. The third is, by 
inspection of dash-dotted lines in Figs. 2(a) and (c), that the shapes of the efficiency curves 
of the modes in wider slits resemble their eigenfunctions |'?/'^"'*^| in Fig. 2(f). The latter 
two features can be attributed to the treatment of MEM which involves a mutual expansion 
between the modes in different layers. 

We should keep in mind that the total field at the layer boundary, Hz{x), can be re- 
spectively achieved by the linear combination of eigenfunctions in the narrower and wider 
slit, and the expansion coefficients depend on the position of the narrower slit, subject to 
boundary conditions. Then, the curves in Figs. 2(a) and (c) can be explained qualitatively. 

We first sec the case where the wave is incident from the narrower slit to wider one, as 
shown by the inset in Fig. 2(a). For reflected waves, the reflection efficiencies are pro- 
portional to the projection 'i/j^"'\x)Hz{x)dx where the integration is along the interface 
between the narrower and wider slits. Since the eigenfunctions in the wider slit or their 
combination, Hz{x), can be seen as a smooth variation within the range of narrower slit, 
the reflection efficiency of the SPP mode in the narrower slit, oc ilj["'°'\x)Hz{x)dx ~ 
J^(na) ip["''^\x)Hz{x)dx, is dominant because its cigcnfunction is also smooth within the slit, 
and that |i?2| oc Jq iJj^"'\x)Hz{x)dx ~ Jyj(na) ip^"'\x)Hz{x)dx is very small because the sec- 
ond mode is an antisymmetric function within the slit. For the transmitted waves, the 
transmission efficiencies |T„| are qualitatively determined by 'il)!il"'^\x)Hz{x)dx, which in- 
cludes /o^V^'^(x)V'i""^(a;)da;, il;^^'\x)i;''r\x)dx , i>^J'\x)i^i"''\x)dx, and so on. We 
have already known that the excitation of the second mode in the narrow slit is very small, 
so that the contribution of the factor ipl^^\x)ilj^"'\x)dx is negligible. The contribution 
of the radiation modes is relatively complicated, but unimportant because what happened 
inside the slit is the key part of the scattering procedure; while the radiation modes local- 
ized in metal are excited to fulflU the boundary condition outside the slit. That is why 
we try to avoid theses modes in the discussion. By several numerical tests, it is sure that 
the radiation modes do have contribution to the transmitted waves but the contribution is 
comparatively small. Therefore, the factor Jq 'iljl^^\x)ilJi^"'\x)dx ^ J^(„a) ip^^''{x)ilj[""'\x)dx 
mainly determines the transmission efficiencies |T„|. As an example, let us see the IT3I curve. 
\Ts\ oc /^(na) i/j2"^\x)xlJi^^\x)dx, where ipi^"'^ is smooth within a narrow region, see, the black 
fine in Fig. 2. When the narrower sht is positioned at the left side with its center being 



10 



Bit X — 0.65 //m, |'03 I has a maximum at this position. Therefore the projection of ipi 
onto ■03"''^ is at a maximum. As the narrower sht moves rightwards, we image that the black 
curve in Fig. 2(e) shifts rightwards. At x = 0.794 //m, iV'a"'*^! is zero. Accordingly, the 
projection of ■j/'j""'' at this position onto as well as IT3I, reaches zero. Between x = 0.65 

and X — 0.794 /im, IT^I should drop from the maximum to zero. We notice that around the 
zero, the phase of T3 changes nearly tt. At the other zero of IV's"'*'*! x — 1.206 //m, 1731 
again reaches zero and its phase changes nearly tt once more. This analysis explains why 
the shape of 1731 is like to iV's""^!- It is the narrow and smooth profile of that causes 

the similarity of the curves between IT3I and iV's"'*''! curves. The IT2I curve in Fig. 2(a) is 
understood in the same way. |Ti| is mainly determined by J^(na) ip^i"^\x)ilj["'"'\x)dx, which 
is a smooth and relatively flat curve due to the smooth variations of both SPP waves. 

We next see the case where a SPP wave is incident from the wider slit to narrower one, 
as shown in the inset of Fig. 2(c), the incident wave being a smooth curve within the range 
of the wider slit width, see the black curve in Fig. 2(f). We again begin with the waves in 
narrower slit. |T„| is proportional to the integral Jq ipl^^"'\x)Hz{x)dx, where Hz{x) is the 
combination of cigenfunctions in the wider slit and considered as a smooth varying curve 
within the range of the narrower slit, so that the transmission efficiency of the SPP mode is 
dominant and much larger than that of the second mode. For the waves in the wilder slit, 
ignoring the contribution of the second mode and radiation modes, the reflection efficiencies 
\Rn\ is mainly determined by f^(na) ^|J^^\x)^jJl^°'\x)dx, leading to the fact that the 
curves in Fig. 2(c) have similar shapes as |T„| curves in Fig. 2(a). 

The transmission efficiency |Ti| in Fig. 2(a) is exactly the same as that in Fig. 2(c), and 
the efficiencies |-Ri|, IT2I and IT3I in Fig. 2(a) have the same behavior as |-Ri|, I-R2I and {R^l 
in Fig. 2(c), respectively, although with different values. The reflection efficiencies \Ri\ and 
|i?2| in Fig. 2(c) are higher than those in Fig. 2(a). This is because Fig. 2(c) represents the 
case that wave incident from a wider slit to a narrower one, which needs to squeeze light 
into a narrower space, so the higher reflection is understandable. 

The variations of the scattering phase shifts for the above two different incident cases 
are plotted in Fig. 2(b) and (d). It is seen that the phases of both Ti in these two figures 
are also the same. We note that at the positions where '0^""^ is zero, the corresponding 
coefficients it!„ and T„ have phase change of tt. 

The results of Type II structure are plotted in Fig. 3. The left walls of the shts are 
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aligned at Xi = 0.6 /iin. The width of the narrower sht is 0.1 /im, but that of the wider one, 
denoted by w, varies from 0.1 to 0.8 yum, as shown in the insets in Figs. 3(a) and (c). The 
most distinct feature is the drastic changes of the excitation efficiencies over a narrow range 
of sht width, as shown in Figs. 3(a) and (c) near w = 0.46 fim. 
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FIG. 3: (color online). Scattering in Type II junction structures. The left walls of the slits are 
aligned at xi = 0.6 /im. 

Q(0) ^ Q(i) ^ Q(2) ^ 0. In all the fi eures, the results of the wider slit 
are plotted by solid lines and those of the narrower slit by dash-dotted lines. The width of the 
narrower slit is 0.1 jiva.. The x-axes are the width of the wider slit, denoted as w. For wave incident 
from the narrower slit to the wider one, (a) excitation efficiency; (b) phase shift. For wave incident 
from the wider slit to the narrower one, (c) excitation efficiency; (d) phase shift. In (e) and (f) 
plotted are the real and imaginary parts of the propagation constant ky^ appearing in Eq. (1) with 
the notation ky \ of the wider slit. 
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When the width of a sht is 0.1 and 0.8 //m, the first three eigenmodes have been plotted 
in Figs. 2(e) and (f), respectively. Now the width w varies. Our calculation shows that the 
modes resemble those in Fig. 2(e) when w < 0.15 /xm and those in Fig. 2(f) otherwise. The 
first three eigenmodes are here labeled by the notations ipi, 1^2 and ip^, and their ky^s are by 
kyi, ky2 and ky^, respectively. 

The SPP mode ipi is obviously a propagation one since kyi has a negligible imaginary 
part, and ips is a decaying one as demonstrated by the large imaginary part of ky^. When 
w < 0.46 /im, ky2 is nearly purely imaginary so that ■02 is an evanescent mode, while when 
w > 0.46 nm, ky2 becomes nearly purely real so that ■02 turns to be a propagation mode. 
A turning point appears at w = 0.46 fim at which ii'2 transforms its propagation property. 
This transformation leads to the drastic changes of the excitation efficiencies, similar to the 
cause of the well-known Wood's anomaly in the grating theory [21]. 

The analysis about the excitation efficiencies in Figs. 3(a) and (c) are in the same way 
as those in the Type I structure. Therefore, some similar conclusions are obtained, such as 
the identity of curves in Figs. 3(a) and (c) and the similarity between the efficiencies 
l-Ril, IT2I and IT3I in Fig. 3(a) and \Ri\, |i?2| and |i?3| in Fig. 3(c), respectively, although 
with different values. 

The drastic changes of excitation efficiencies at the turning point have to be explained 
from an energy perspective. As an example, let us see the case where the wave is incident 
from the narrower slit to wider one, as shown by the inset in Fig. 3(a). At the start point 
w — 0.1 //m, the two slits are the same, so that |Ti| = 1 and all other excitation efficiencies 
are zero. Close to the turning point, |Ti| reaches the minimum, and and IT2I reach the 
maximum. Since the second mode in the wider slit 1^2 now is an evanescent mode, the energy 
is mainly stored in the reflected SPP wave. Once ip2 becomes a propagation mode, it must 
gain a large portion of energy from the reflection due to its large excitation efficiency, and 
leads to the rapid drop of \Ri\ as shown in Fig. 3(a). At the same time, the dropping 
further causes a redistribution of excitation efficiencies by the boundary continuum condi- 
tion. Therefore, what behind the drastic changes of excitation efficiencies is a redistribution 
of energy between evanescent modes and propagation modes. The explanation of the energy 
redistribution is also suitable for the case where the wave is incident from the wider slit to 
narrower one. 

It is seen from Figs. 3(b) and (d) that the phases of both Ti in these two figures are also 
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the same. At the turning point, Ti changes its phase by tt. The change of tt in phase at the 
turning point also occurs for Ts in Figs. 3(b) and R3 in Figs. 3(d). 

So far, the scattering mechanisms of the two types of structures are investigated. Al- 
though the investigation here is restricted to the incident wave with wavelength Aq = 1 /xm 
only, the analysis above is also applicable to the infrared and visible spectrum. Furthermore, 
the analysis is important in practical application. For example, one can excite/suppress spe- 
cific modes to control the field distribution inside a sht, or design high efficient reflector, by 
changing the position or width of the slit. 

B. The multi-mode mult i- reflect ion model 

Having had the knowledge of the scattering of the interface in a single junction structure, 
we are ready in this subsection to discuss the transmission in a step-modulated slit which 
can be regarded as the combination of two junction structures. In order to reveal the 
transmission clearly, we present here an analysis of multi-mode multi-reflection model that 
combines wave and ray optics. 

Figure 4 is the sketch of the multi-mode multi-reflection model. In a step-modulated 
subwavelength metal slit, there are two interfaces at Q^^^ and Q*-^^ Let us discuss the wave 
reflection and transmission in the slit. When the incident SPP launched from Q^^^ in Layer 1 
impinges the interface between Layers 1 and 2, Q^^\ it generates the reflection wave in Layer 
1 and transmission wave in Layer 2. The latter continues going upwards, and when reaching 
other interface Q^'^\ yields reflection and transmission waves again. Obviously, there occurs 
multi-reflection in Layer 2, shown by Fig. 4(a). 

Because the wave in each layer is the linear combination of the eigenmodes of the layer, 
each ray in Fig. 4(a) can in fact be expanded by eigenmodes, except the primary incident 
light. For example, at the point A, the reflected wave contains all possible modes in Layer 
1 and the transmitted wave contains the modes in Layer 2. That is to say, the scattering 
excites all the modes in both layers. When all the possible modes in Layer 2 reach point 
B, each mode again excites all possible eigenmodes in reflected wave in Layer 2 and in 
transmitted wave in Layer 3. The phenomenon is termed as multi-mode excitation. To 
show the phenomenon explicitly, we draw in Fig. 4(b) the multi-mode excitation at point C. 
Suppose that the waves in Layers 1 and 2 are expanded by three eigenmodes, respectively. 
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FIG. 4: (color online). Sketch of the multi-mode multi-reflection model, (a) multi-reflection 
between interfaces and (b) multi-mode excitations at each point. 

Then when the three modes in Layer 2 are incident to point C, as shown in Fig. 4(b), the 
first mode yields the reflected and transmitted waves, both containing three eigenmodes 
in respective layer, i.e., the incident black line excites the black, red and blue lines in the 
transmitted waves in Layer 1 and reflected waves in Layer 2, respectively. In the same way, 
the incident red line also excites the black, red and blue lines in the transmitted waves 
in Layers 1 and reflected waves in Layer 2, respectively, and so does the incident blue 
line. Therefore, the total reflected wave at point C includes three eigenmodes in Layer 2, 
each being in turn the superposition of the reflections from the three incident eigenmodes. 
Similarly, the total transmitted wave at point C includes three eigenmodes in Layer 1, each 
being in turn the superposition of the transmissions from the three incident eigenmodes. 

In summary, the total transmission and reflection coefficients in Layer 3 and Layer 1 in 
Fig. 4(a) are obtained by summing up all the single-scattered coefficients, respectively. In 
addition, it is worth mentioning that the multi-mode multi-reflection model is a generalized 
form of the single-mode multi-reflection model which occurs in a F-P cavity. The former 
will be simplifled to be the latter if only one mode can be excited. 

The physical explanation of the multi-mode multi-reflection process is named as model 
analysis. In order to testify this analysis, numerical calculation based on this physical picture 
is carried out and the results are compared with MEM. In Fig. 5 plotted are the transmission 
efficiencies as a function of the length of Layer 2 g*^^^ when the incident wave is SPP mode. 
The solid lines in Fig. 5(a) and (b) are the results of MEM, which surely comprise the 
contributions from all possible eigenmodes. The symbols are the results from the model 
analysis. In a slit with width w'^'^^ = 0.3 /im, only the first mode, i.e., the SPP mode, can 
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propagate and all the other modes are evanescent. Thus, when g^^) is sufficiently long, the 
higher modes attenuate to a negligible value, and the transmission can be well described by 
a multi-reflection of only the SPP mode. 
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FIG. 5: The SPP transmission vs. g^^^ calculated by MEM and model analysis for slit structure with 
[q(°),Q(^),(5(2) = -1,0, qC^) i_im. The solid lines are the results of MEM. The circles, crosses 
and plus signs are the results including contributions from the first one, two and three modes, 



respectively, from the model analysis, (a) Slits align to left at = 0.85 /xm, 
[0.1,0.3,0.1] nm; (b) The symmetric case of the slit structure in (a). 
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In Fig. 5(a) it is seen that the results of the model analysis including the first two modes 
are accurately the same as the hne from MEM. When g^^^ > 0.8 iim, the circles and crosses 
are identical, indicating that the contribution from the second mode is negligible. While for 
g*^^) < 0.8 fxm, crosses deviate from circles, indicating that the second mode should not be 
omitted since it does not fade out within this distance range. In Fig. 5(a), the crosses end 
at g*-^-* = 0.09 /im, because below this distance the multi-reflection of the flrst two modes 
diverges. What is the reason of the divergence? Firstly, the divergence is not caused by the 
propagation mode since the multi-reflection of the SPP mode always converges, as shown by 
the circles in Fig. 5(a). Secondly, it is neither caused by the exponentially increasing term 
which originates from improperly handling the evanescent waves [20] for it occurs only at 
short distances. Actually, the divergence arises from the coupling between the eigenmodes 
containing the evanescent modes. With the contribution of the evanescent wave, as shown 
in Fig. 4(b), the superposition will result in a larger transmission and reflection coefficients 
after each scattering if the second mode does not decay to a certain value. Thus, there exists 
a critical distance above which the multi-mode multi-reflection analysis is applicable. For 
the structure given in Fig. 5(a), it is g^^^ = 0.09 iim. 
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To verify the above conclusion about the divergence of the evanescent mode, we suppress 
the antisymmetric second mode by reforming the sht structure to a symmetric one, as 
shown in Fig. 5(b). Then the circles and crosses are identical at any distance. Let us see 
the contribution from the third mode. The plus signs including contributions from the first 
three modes are in good agreement with the results of MEM as q'^'^^ > 0.012 fim. The crosses 
and circles are identical when q^'^^ is above 0.3 /xm, but it is not so when g^^^ is below 0.3 
/im. That is to say, if the distance is less than 0.3 /im, the evanescent third mode is not 
negligible. This time the critical height for the third mode is q^"^^ — 0.012 //m, which is much 
smaller than that of the second mode in the structure shown in Fig. 5(a). The reason is 
that the decay of the third mode is faster than that of the second mode, for the propagation 
constant ky of the former has a larger imaginary part than the latter, as shown in Fig. 3(f). 
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FIG. 6: The SPP transmission vs. g^^^ calculated by MEM and model analysis for slit structure 



with 



-1,0, g*-^^ fim. The solid lines are the results of MEM. The circles, 
crosses, plus signs, up-triangles and down-triangles are the results including contributions from the 
first one to five modes, respectively, from the model analysis, (a) Slits align to left at xf^ = 0.6 
= [0.1,0.8,0.1] fim; (b) The same structure but the narrower slits are shifted 
0.744 jim, a position which totally suppresses the excitation of the 3rd mode. 



^(1) (2)_,„(3) 



to be X 



(1) _ ^(3) _ 



see Fig. 2(a). 

Next, we investigate the coupling between propagation modes. In Fig. 5(a), only the 
first mode, the SPP mode, is the propagation one in Layer 2 with width w^"^^ = 0.3 fim. 
When the width is enlarged, the second mode can also become propagating. In Fig. 6(a) 
plotted are the transmission efficiencies in the same structures as in Fig. 5(a) except that 
the width of Layer 2 is extended to be 0.8 //m. Under this width, the second mode is indeed 
propagating. It is seen that even q^^^ gets to zero, the result containing the contributions 
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from the first two modes is not divergent, and agrees with the MEM curve very well when 
g'(^) > 1.2 fj,m, which confirm the statement that the propagation modes do not cause the 
divergence. When g*^^) is below 1.2 /xm, the contribution from the third mode has to be 
added in order to achieve precise results. However, as the cost of preciseness, the divergence 
appears below q^"^^ = 0.181 fim. 

Similar to the treatment in Fig. 5(b) where the second mode in Layer 2 is removed by 
structural change, it is also possible to suppress the third mode excited in Layer 2. The way 
to implement the suppression is to shift the center of the narrower slits to x — 0.794 /im, as 
shown in the inset of Fig. 6(b). At this position, the excitation efficiency of the third mode 
is nearly zero, see. Fig. 2(a). The transmission results are plotted in Fig. 6(b). It is seen 
from the figure that up to the first five cigenmodes have to be included in the model analysis 
in order to meet the MEM curve. The divergence in this case is caused by the fourth mode, 
and the corresponding critical width is g^^^ = 0.084 //m. 

In summary, the multi-mode multi-refiection model provides intuitive and precise de- 
scription about the transmission inside a step-modulated subwavelength metal slit when the 
height of modulated layer (Layer 2) is above a critical height, while fails below it because of 
the coupling between propagation modes and evanescent modes. 

C. Comparison of different methods 

In this subsection, MEM and other three methods, FDTD, SMM, and ICIM, are dis- 
cussed, and the calculated results of FDTD and ICIM are compared to the MEM results. 
The preciseness of these methods is investigated and some useful conclusions are obtained. 

Before presenting the numerical results, we would like to make a brief discussion about these 
four methods. 

FDTD, as a commonly used simulation method in optics, is to calculate field quantities 
directly from the Maxwell's equations by difference method. In principle, this method and 
MEM both can provide accurate and reliable results. Here we would like to point out 
their three discrepancies. Firstly, the way they solve the Maxwell's equations is different: 
FDTD uses finite difference method to evolve fields in space and time domains, while the 
MEM establishes and solves the coupled equations in frequency domain by the method of 
moments. Secondly, the way they handle outmost boundaries is different: FDTD makes 
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use of, for the outmost boundaries of a system, perfectly matched layers which can totally 
absorb waves without reflecting them back, while MEM confines the structure with two 
perfectly conducting walls such as in this paper. The feasibility of the latter is due to the 
fast attenuation of light (infrared, visible spectrum) in a metal. If the confined width is large 
enough, the effects brought by the two perfectly conducting walls are negligible, as shown 
in Fig 7(a) below. Although the perfectly matched layer technique can be introduced to 
MEM [22], it dramatically complicates the modal analysis. Thirdly, the way they converge 
is different: the convergence of FDTD depends on the size of the Yee cell used in simulation, 
while that of MEM on confined width L and truncation number N . 

SMM [7-9] and ICIM [13,14] are other two frequently used methods which show following 
three features. Firstly, according to the two methods, the modulated region. Layer 2, would 
be divided into a central scattering region and a stub (as shown in the Fig. 2 in Ref. [7] 
and Fig. 4 in Ref. [13]), and it would assume that the SPP mode multi-reflection occurred 
in the stub (although Refs. [13] and [14] did not mention this point, it could be recognized 
from the transmission equations, Eq. (4) in Ref. [13] and Eq. (8) in Ref. [14]). Secondly, 
both of them took the one mode approximation, which meant that only the SPP modes 
existed in the stub and slits. Thirdly, the phase shifts caused by scattering in the central 
scattering region could not be calculated properly, so that were ignored by means of the 
quasi-statistic approximation [23,13,14] (although Refs. [7-9] did not mention the quasi- 
statistic approximation, it was easily seen by the procedure of obtaining scattering matrix 
given in Ref. [8]). 

In the following, the numerical comparison between these methods is performed. The 
convergence comparison of MEM and FDTD in a Type I structure is presented in Fig. 7. 

Figure 7(a) shows the calculated results of MEM with confined width being L = 2, 4 jim 
and truncation number being N — 400, 800, respectively. The results by squares, crosses 
and plus signs in Fig. 7(a) are identical, showing that the boundary effect imposed by the 
perfectly conducting walls can be completely ignored in such confined widths. As already 
mentioned at the beginning of Sec. 3, the parameters L = 2 /im and N = 800 ensure 
that all the calculated results have at least four significant digits. A FDTD simulated 
transmission curve with the cell size being 2.5 x 2.5 nm^ is also plotted in Fig. 7(a) for 
comparison. Obviously, the FDTD curve is close to the MEM ones but not coincide. In 
order to investigate the convergence of FDTD, three simulated transmission curves with 
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FIG. 7: (color online). The SPP transmission and reflection in a Type I structure with parameter 
= 0.6 /im, U(i),u;(2),u;(3)l = [0.1,0.8,0.1] //m, [q(°), Q(^), Q^^)! = [-1,0,0.1] ^m, see the 
inset of (a), (a) Convergence test for MEM; (b) convergence test for FDTD relative to MEM 
where x-axis is a part of (a) for xj^^ G [0.7, 0.8] 

different cell sizes are plotted in Fig. 7(b) for the same structure as in Fig. 7(a) but with 
horizontal abscissa being x^i^ G [0.7, 0.8] /xm to highlight two absorption peaks. A MEM 
curve also plotted in the figure for comparison. For large cell size as 5 x 5 nm^, only one 
vague dip, instead of two, is observed. The dip will gradually separate into two and approach 
to the MEM curve as the cell size decreases. However, even for 1.25 x 1.25 nm^, namely, 
1/800 of the incident wavelength or 1/80 of the stub width (height of Layer 2), the deviation 
of the results between FDTD and MEM is still observable, which means that FDTD has 
a relatively slow convergence. That is why we do not use FDTD to verify the calculated 
results of MEM in this paper. 

The transmission of the structure considered in Fig. 7 was also investigated by SMM 
[9]. This method actually utilizes some results of FDTD to obtain scatting matrix elements 
and loses some phase information by quasi-statistic approximation, so that its final results 
could not be better than that of FDTD. This can be recognized by the comparison between 
the results of FDTD and SMM given in Refs. [8] and [9] . Therefore, it is not necessary to 
discuss the preciseness of SMM here because it depends on the simulation results of FDTD 
which has already shown in Fig. 7(b). Besides, since the SMM and ICIM have the similar 
transmission model, their calculation errors ought to have the same order of magnitude. The 
calculation error of ICIM is investigated in the following. 

The SPP transmission calculated by ICIM and MEM in a Type II structure are plotted 
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in Figs. 8 and 9. This kind of structure was also study in Refs. [13] and [14]. 




FIG. 8: Comparison of the SPP transmission by MEM and ICIM of a Type II structure under 
the variation of w^"^^ and q^'^\ The left sides of the slits in ah layers are aligned to xf^ = 0.6 
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by MEM; (b) \tI'^i^\ by ICIM; (c) absolute value of the difference between MEM and ICIM, 
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for g(^) < 0.46 /xm. 



In Fig. 8 are given the SPP transmission as a function of the height and width of Layer 
2 under a fixed incident wavelength, Aq = 1 /im. In Fig. 8(a), the calculated results of 
MEM can be approximately divided into three areas A, B, and C by dotted lines. In Area 
A, a regular oscillating pattern is observed because when w^"^^ < 0.4 /im only one SPP mode 
propagates in Layer 2 that forms the FP-hke oscillation [5,6]. While in Areas B and C 
where the higher modes also contribute to transmission, the transmission pattern becomes 
complicated. Intuitively, there are two ways for higher modes to transport energy. One is 
in a way of a propagation mode, which is appropriate for w^^^ > 0.46 fim because the 2nd 
mode in Layer 2 becomes propagating. The other is in a way of an evanescent mode, which 
is appropriate for 0.4 < w'^'^^ < 0.46 /xm and Area C because the 2nd mode in Layer 2 does 
not attenuate to a negligible value in such a modulated slit and brings energy through Layer 
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2. 

The results from ICIM shown in Fig. 8(b) can be discussed according to if g^^^ is larger 
or less than 0.46 /xm. When q^'^'> > 0.46 fim, more than one mode are allowed to propagate 
in the stub, while the ICIM assumes only the SPP mode, which raises the great difference 
relative to the MEM results, leading to totally different patterns between Figs. 8(a) and 
(b). When < o.46 /im, the difference \Tf^E'^\ - \tI^'im^ jg plotted in Fig. 8(c). In 
this region, the difference is mainly caused by the neglect of the phase shifts in the central 
scattering region, which indicates that these phase shifts have to be taken into account in 
calculation. In Fig. 8(c), it is seen that when the length of Layer 2 g^^^ < 0.03 //m, less than 
1/30 of the incident wavelength, the difference becomes larger due to the energy transported 
by evanescent modes. This demonstrates that a very narrow region cannot guarantee precise 
results. Thus, the application of the quasi-statistic approximation in a modulated metal slit 
requires an optimum geometry in order to provide relatively accurate results: the incident 
wavelength is nearly 10 times larger than stub width. 

In order to test the optimum geometry, the slit width of Layer 2 and the incident wave- 
length Ao are considered as variables and the length of Layer 2 is set as g*^^^ = 0.1 /xm, 
the same as the slit widths in Layers 1 and 3. The MEM results in Fig. 9(a) exhibit a few 
straight black bands indicating the inverse-proportional relation between frequency and stub 
width which was mentioned in Refs. [8] and [9]. For ICIM, similar results are obtained in 
Fig. 9(b). The difference of these two figures is plotted in Fig. 9(c). Clearly, the calculation 
error of ICIM in this case is lower than that in Fig. 8(c). Although small deviation takes 
place near the absorption peaks, ICIM successfully predicts the positions of the peaks. All 
these numerical results confirm that the ICIM can provide relatively accurate results when 
the incident wavelength is around 10 times of the stub width. 

Anyway, each method has its own advantages in certain aspect. For example, FDTD can 
provide visualized transmission process in time domain and ICIM is indubitably the fastest 
method in calculation. Here we just emphasize that these methods should be used with 
caution in considering the convergence and preciseness. 
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FIG. 9: Comparison of the SPP transmission by MEM and ICIM of a Type II structure under 
the variation of uP''^ and wavelength Aq. The left sides of the slits in all layers are aligned to 
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IV. CONCLUSION 

In this paper, the MEM developed in Ref. [6] is improved to be a more practical and 
efficient one for handling the scattering/ transmission. Using this method, the scattering in 
a juncture structure and the transmission inside a step-modulated slit are investigated. 

Firstly, the scattering in a juncture structure is studied for two types of structural changes. 
For the Type I change where the widths of the two shts are fixed and the position of the 
narrower one can be anywhere within the wider one, the excitation efficiencies of the modes 
in the wider slit resemble their eigcnfunctions respectively, while in the narrower slit the 
excitation efficiency of the SPP mode is dominant and much larger than that of the second 
mode. For the Type II change where the left walls of the shts are ahgned and one sht 
becomes wider gradually, a wood-anomaly-like drastic change of excitation efficiencies is 
observed when the propagation property of one mode begins to transform. 
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Then, the transmission inside a step-modulated sht is studied. Besides the MEM calcu- 
lation, we present explicitly a multi-mode multi-reflection model to reveal the transmission 
process. The multi-mode excitation and the superposition procedure of the scatterings from 
all possible modes are the key parts of the model, which represent the interference and en- 
ergy transfer happened at layer boundaries. However, there exists a critical height of the 
modulated layer for applying the model due to the coupling between propagation modes 
and evanescent modes. Above the critical height, the model can provide the same result as 
MEM. 

In addition, some commonly used methods are compared with MEM. Useful conclusions 
about these methods are obtained: for a subwavclcngth metal slit, MEM is a versatile and 
fast method that can provide accurate results; FDTD has a relatively slow convergence and 
need very small Yee cell to ensure the accuracy of the simulation; ICIM incorporating the 
one mode and quasi-statistic approximations provides relatively accurate results when the 
incident wavelength is around 10 times larger than the stub width. 
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